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Most of the literature on quantum vortices predicting various states of vortex matter in three 
dimensions at finite temperatures in quantum fluids is based on an assumption of an extended and 
homogeneous system. It is well known not to be the case in actual Bose-Einstein condensates in 
traps which are finite systems with nonuniform density. This raises the question to what extent 
one can speak of different aggregate states of vortex matter (vortex lattices, liquids and tensionless 
vortex tangle) in these system. To address this point, in the present work we focus on the finite- 
size, boundaries and density inhomogeneity effects on thermal vortex matter in a Bose-Einstein 
condensate. To this end we perform Monte Carlo simulations on a model system describing trapped 
Bose-Einstein condensates. Throughout the paper, we draw on analogies with results for vortex 
matter obtained for extended systems. We also consider, for comparison, the cylindrical confinement 
geometry with uniform density profile from the center out to the edge of the trap. The trapping 
potential is taken to be generically of an anharmonic form in such a way as to interpolate between 
a harmonic trap and a cylindrical confinement geometry. It also allows for a dip in the density 
profile at the center. We find distinct thermal equilibrium properties of the vortex system as the 
qualitative characteristics of the trapping potential is varied. For a uniform cylindrical confinement 
geometry, we find a vortex lattice at the center of the trap as well as ring-like thermal fiuctuations 
of the vortex system as the trap edge is approached. For a harmonic trap, we find a distinct region 
at the edge of the trap where the vortex lines appear to have lost their line tension. Due to the 
steep density gradient, this crosses directly over to a vortex-line lattice at the center of the trap at 
low temperatures. At higher temperatures, an intermediate tensionful vortex liquid may exist. For 
an anharmonic trap where the ground state condensate density features a local minimum at the 
center of the trap, we find a corresponding region which appears to feature a tensionless vortex- 
line liquid phase. This work suggests that finiteness and intrinsic inhomogeneity of the system not 
withstanding, one nonetheless can approximately invoke the notion of distinct aggregate states of 
vortex matter realized at certain length scales. This might be helpful, in particular in search of 
possible new states of vortex matter in Bose-Einstein condensates with multiple components and 
different symmetries. 



PACS numbers: 03.75.Hh,03.75.Kk,67.40.Vs 

I. INTRODUCTION 

The physics of trapped gases on the one hand, and the 
physics of superconductors and superfluids on the other, 
may be conceptually linked by rotating Bose-Einstein 
condensates (BEC) in magnetic traps^, or pair conden- 
sates of ultracold Fermi gases^^^^^. Superconductivity, su- 
perfluidity, and Bose-Einstein condensation are all hall- 
marks of quantum fluids governed by the onset of long- 
range phase coherence in a macroscopic matter wave. 
This phase-coherence is the matter- wave analog of a cor- 
responding well-known phenomenona in electromagnetic 
waves, namely phase-coherence in such waves established 
by stimulated emission of radiation. One distinguishing 
feature of such quantum fluids is that the macroscopic 
matter wave function is complex ip{r) = |?/^(r)|e*^^^\ and 
with a phase 9{r) G [0, 27r), i.e. this phase is defined with 
a compact support. This has far-reaching ramifications 
for the physics in the sense that the order parameter of 
the system supports stable topological defects in the form 
of point- vortices in two dimensions and vortex-loops and 



vortex-lines in three dimensions. Phase transitions from 
superfluids to normal fluids, are entirely governed by such 
topological defects. 

In two dimensions, this is manifested in the well known 
Kosterlitz-Thouless phase transition of unbinding of a 
vortex-antivortex pair from a low-temperature superfluid 
with only tightly bound vortex-antivortex pairs to a nor- 
mal fluid in a process where the most weakly bound pairs 
are dissociated^. It is also equivalent to a phase transi- 
tion from a dielectric to a metal in the two-dimensional 
Coulomb-gas with overall charge-neutrality^. In three 
dimensions, one has a phase-transition from an ordered 
low temperature phase with at most a few small closed 
vortex loops present, to a normal fluid where closed vor- 
tex loops have proliferated throughout the system into a 
tangle in such a way as to make it possible to connect 
opposite sides of a macroscopic system with a connected 
vortex pathij^i2ii^5^ii^. This was originally proposed by 
Onsager already in 1949 in qualitative terms as a way of 
explaining the A-transition in superfluid ^Hei^. In the 
context of the present paper, we stress that while Refs. 
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I ll3l dealt with the issue for zero rotation for the super- 
fluid or zero magnetic fleld for the superconductor, the 
picture was extended to flnite rotation in superfluids and 
flnite magnetic fleld in extreme type-II superconductors 
in Refs. 8,9,10. It is the latter situation that is relevant 
to rotating Bose-Einstein condensates. 

The connectivity of the thermally excited vortex tan- 
gle changes abruptly at the critical temperature of the 
system. Such a phase transition in three dimensions can- 
not be cast into the framework of a Kosterlitz-Thouless 
phase transition for the simple reason that vortex loops in 
three dimensions have a property that vortex- ant ivortex 
pairs in two dimension do not have, namely shape. This 
contributes signiflcantly to the conflgurational entropy 
of the system. The fact that vortex loops can have ex- 
tremely complicated geometric shapes and will form a 
fractal structure at long length scales at the critical point, 
is crucial in order for them to be able to proliferate. It 
also means that the vortex loops cannot be regarded as 
renormalized vortex rings with a 'doughnot hole' in the 
middle. They are instead fractal objects with fractal di- 
mension Dh ^ 2^. Such a fractal dimension is consid- 
erably larger than what it would have been for ring-like 
objects, namely Dh = 1. The entropy production as- 
sociated with this proliferation of topological defects is 
accompanied by a loss of a generalized stiffness, in this 
case the superfluid density or phase stiffness of the sys- 
tem. In this sense, the above scenarios both in 2D and 
3D fall nicely within a general deflnition of a phase tran- 
sition, namely that a phase transition occurs at a point 
where some generalized stiffness is lost as a result of a 
spontaneous proliferation of stable topological defects of 
the complex scalar matter flelds in the system^^. 

Over the last decade, remarkable progress has been 
made in achieving Bose-Einstein condensates in gases of 
ultra-cold atoms in a magnetic trap^^^^*"'^^. Such con- 
densates are now being routinely manipulated in a large 
variety of ways, and may for instance be spun up to 
produce vortex lattices of a condensate in a magnetic 
trapi^i2£. One may also envisage low-dimensional vor- 
tex structures^^^^^. There are even cases where other 
aggregate states of vortex matter are known to exist in 
quite different condensed matter systems, such as two- 
component superconductors with individually conserved 
condensates^^. It is the purpose of this paper to study the 
thermal excitations of vortex-lines and vortex-loops in 
Bose-Einstein condensates which are conflned to a cylin- 
drical geometry by a trapping potential. This trapping 
potential generally increases from the center of the trap 
towards the edge of the condensate, although this vari- 
ation may be non-monotonic. In essence, it acts as a 
spatially dependent chemical potential for the bosons in 
the system, thus affecting the overall condensate density. 
The density is typically highest at the center of the trap 
and vanishes towards the edge of the trap, although more 
complicated proflles may easily be envisaged, and will 
in fact be considered in this paper. The way in which 
the overall density varies is directly determined by the 



trapping potential. Thus, such systems are inherently 
nonuniform and therefore it is important speciflcally to 
study the effect of spatial density variations. 

This paper is organized as follows. In Section [III we 
motivate and introduce the model on which we will per- 
form Monte Carlo simulations in this paper. In Section 
mil we describe the Monte Carlo simulations we will per- 
form. In Section IIVI we present results for a uniform 
system, as a benchmark, and we also consider such a 
system conflned in a cylindrical geometry. The former 
of these results are applicable to fluctuating vortex mat- 
ter in extreme type-II superconductors, where we can 
neglect the fluctuations in the gauge fleld. The latter 
results apply to rotating helium in a container. In Sec- 
tion [Vl we present our results for two different types of 
traps, namely the harmonic cylindrical trap, and the an- 
harmonic cylindrical trap. The quantities we focus on 
are the helicity modulus and structure functions of the 
vortex matter in the systems in the various parts of the 
conflnement and traps, i.e. at various distances from the 
center. In Section I VII we present our conclusions. This 
work is a follow-up of the letter Ref. i25,. 



II. GENERAL MODEL 

We will use Monte Carlo (MC) computations to search 
for the thermal equilibrium vortex states of a rotating 
Bose-Einstein condensate in three dimensions. A natu- 
ral starting point is the Gross-Pitaevskii free energy func- 
tional in a rotating frame, given by the energy functional 



dV 



^ 2M 



-V{r)-nL,)^^-g\^\^ 



(1) 



Here, Lz = \h{ydx — xdy) and g = Aiiash^ /M ^ M is the 
mass of the atoms in the condensate, and is the 5- wave 
scattering le ngth. The condensate wave function is given 
by il){r) = y^^p(r)e'^^^^\ where p{r) is the local density of 
the bosonic matter in the trap. Long-range ordering in 
the phase variable 0{r) signals the onset of superfluidity 
or Bose-Einstein condensation. It will be important for 
our later discussion that, consequently and conversely, 
the destruction of the Bose-Einstein condensate proceeds 
via phase-disordering of the system through large phase 
fluctuations V^. Here, "large" means phase fluctuations 
that give rise to vortices due to the compact nature of 
the phase variable 6 G [0,27r), i.e. phase fluctuations 
that have the property V x (V^) = 27rn, where n is 
an integer- valued vector. These are the transverse phase 
fluctuations, as opposed to the longitudinal ones, which 
correspond to spin- waves in an XY ferromagnet. 

The potential V{r) is a trapping potential which con- 
flnes the Bose-Einstein condensate in a flnite region in 
space. Moreover, it is also seen to appear as a spatially 
varying chemical potential for the condensate density and 
as such will set the overall density proflle of the con- 
densate. Increasing V{r) will suppress the condensate 
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density and vice versa. Furthermore, Q is the angular 
frequency associated with the rotation of the condensate 
when it is spun up, and Lz is the corresponding total an- 
gular momentum operator. Finally, the quart ic term is 
a repulsive contact interaction between the bosons of the 
condensate which will render the spectrum of the theory 
bounded from below. Such an energy functional is ap- 
plicable as a coarse grained description of an uncharged 
phase coherent condensate^. 

If we formally introduce a vector potential A = 
{M/h){fl X r) with ft = l^z, the energy functional in 
Eq. ([1]) can be written 



dV 



2M 



^*(V-iA)2^ 



{v{r)--Mn\l)m' + -g\^\ 



(2) 



where 



y'^. This is formally similar to the 



Ginzburg-Landau free energy of a superconductor, apart 
from the position-dependent term involving the trapping 
potential V{r) and rotation Q. However, in a supercon- 
ductor the vector potential A has dynamics and is related 
to the magnetic field B by B = V x A, and thus we would 
normally include a term in the free energy 

density. For extreme type-II superconductors where the 
typical penetration depth Al of a static magnetic field 
is much larger than the coherence length ^, fiuctuations 
in the magnetic field can in many cases be ignored and 
this term in the energy can be dropped. The similarity 
of the superconductor Ginzburg-Landau theory to that 
of a rotating Bose-Einstein condensate is then striking. 

So far, the energy functional of Eq. ([2]) is a contin- 
uum theory, but for the purposes of carrying out com- 
puter computations, it is more convenient to discretize 
space into a lattice and to rescale the wave function 
to avoid prefactors. That is, we let the wave function 
ip{r) [s/M /fi)il)i^ so that it is only defined on vertices 
i = 1, . . . , L^, separated by a lattice constant a. We also 
replace the gradient term with a gauge invariant lattice 
difference. 



■iAi 



(3) 



i + a/i is the lattice site situated next to site i in direction 
/i, and the gauge field Ai^ here lives on the links of the 
lattice and it is given by the line integral 



(4) 



The continuum theory is recovered if we let a ^ 0. 

For high-Tc superconductors, where Xl ^ ^, it is a 
well established approximation only to consider fiuctua- 
tions in the phase of the order parameter. This is fre- 
quently referred to as the London approximation^^, in 
which we simply assume a condensate of Cooper pairs 
to exist by having a finite and fixed |?/^^| = |?/^|, since in 



the end, it is not the depletion of the number of Cooper 
pairs that is responsible for destroying superconductiv- 
ity, but rather the proliferation of vortex loops originat- 
ing with violent transverse phase fiuctuations in the su- 
perconducting order parameter^' High-Tc supercon- 
ductors are extreme type-H superconductors, which in 
a certain sense means that the (renormalized) charge 
of the condensate is small^*^. This suppresses gauge 
field fiuctuations and an external magnetic field there- 
fore only acts as a frustration on the system via mini- 
mal cupling to a fixed external vector potential, just as 
in the above Eq. ([2]). Superfiuids and Bose-Einstein 
condensates have zero charge and may in this sense be 
viewed as the ultimate extreme type-H superconductors 
where all vestiges of the fiuctuating gauge field in the 
problem have vanished. In these systems we therefore 
expect that the phase only approximation is essentially 
exact. Increasing the temperature, large fiuctuations in 
the phase makes the condensate incoherent and non- 
superconducting before \ip\ vanishes, see for instance Fig. 
1 of Ref. [TO. In fermionic pair condensates^^i^, it has 
been strikingly demonstrated that one may have pairing 
without superfiuidity^^'— providing further confir- 
mation of the above ideas in a completely different setting 
than extreme type-II superconductors or ^Hei^. 
The right hand side of Eq. ([3]) can be rewritten 



|tA|2[2-2cos(A^e, 



Ai^.)], (5) 



using the lattice difference operator A^6>i = ^i+a/i — ^i- 
A slight generalization is instead to replace V^i+a/i and 
tjji in Eq. ([3|) with their average as this will allow for a 
non-uniform condensate density^ which is what we want 
to study here. 

Since we will only include phase fiuctuations in the 
Monte Carlo computations, any terms not containing 
will represent mere constant shifts in the total energy 
and we will consequently drop them. We thus arrive at 
a simple, effective energy 



E ■ 



^ ^i/i cos(A^6>i - A^^). 



(6) 



The sum is over all positions and directions x, and 
except for a factor M/h^^ the position dependent cou- 
pling Pi^ is nothing but the condensate density at the link 
from site i in /i-direction. If this factor is set to unity, Eq. 
([6j) reduces to the well known uniformly frustrated 3D 
XY model, used for modelling the melting of the vortex- 
line lattice in uniform bosonic condensates and extreme 



type-II superconductors (see e.g. Refs. I34ll35ll36ll37ll38h . 
The local vorticity can be calculated from a phase config- 
uration by summing the gauge- invariant phase difference 
around each plaquette in the numerical grid. 



(7) 



Here, nj is the number of vortices penetrating a plaquette 
□j in the positive direction. The spatial (radial) varia- 
tion of Pi^ refiects directly the spatial variation in the 
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ground state of the system of the quantity |'?/^(r)p due to 
the spatial variation of the effective chemical potential 
V{r) - ^Mn'^rl_ appearing in Eq. ([2j). 

We mention in passing that longitudinal phase fluc- 
tuations are innocuous in 3D, and hence need not be 
considered for the purposes of studying phase transitions 
in the system. Such phase-fluctuations are incapable of 
destroying the superfluid density in 3D. In 2D they suf- 
fice to render the system critical at any finite tempera- 
ture below the Kosterlitz-Thouless transition, meaning 
that phase-correlations ^(r — r') = ^e*^^'^^e~*^^'^ ^) ex- 
hibit a quasi- long range power-law decay G(r) ~ 
with a temperature-dependent exponent rj. In neither 
case are longitudinal phase fluctuations capable of driv- 
ing the system through a phase transition and into a 
phase with short-range, exponentially decaying phase- 
correlations G{r) ~ e"^/*^, where ^ is the correlation 
length. 

It should be noted that the model Eq. (j6|) does not 
apply to the Lowest Landau Level regime (LLL) of an 
atomic Bose-Einstein condensate. The latter can be 
identified through the ratio of the interaction energy scale 
to the level spacing of the transverse harmonic confine- 
ment in an axially symmetric trap. 



The result is a Markov chain of configurations that can 
be used to estimate the partition function 



MhuJi_ ' 



(8) 



where n is the particle number density, is the s-wave 
scattering length, M is the particle mass and u;_\_ is the 
trap frequency. The LLL approximation is considered 
to be valid when A <C 1^^^^^. However, we expect Eq. 
([6]) to be adequate under the conditions of those ex- 
periments which meet the followi ng naiv e requirement. 
The intervortex distance 2ro = 2^h/MVt should be sub- 
stantially larger than the healing (or coherence) length 
^ = h/Mu;±R±^ where R± is the radius of the system 
perpendicular to the axis of rotation^^ i^^i^^i^^i^^i^^ . 



III. MONTE CARLO COMPUTATIONS 

The Monte Carlo computations are performed with the 
standard Metropolis-Hastings algorithm^^^^^, and the up- 
dates are always local. Initially, all phases are chosen 
equal to zero, but any other configuration would suffice, 
provided that the system is allowed to thermalize for a 
sufficiently long time. We then proceed systematically 
through all lattice sites one by one and propose trial up- 
dates of the local phases Oi 0i-\-50i^ where SOi is drawn 
from a uniform distribution [— 7r,7r). Each trial update 
is accepted with a probability p determined from the en- 
ergy difference of the configurations before and after 
the update. 



{0} 



-B/T 



(10) 



where the sum is over all possible sets of the phase. An 
obvious effect of generating phase configurations accord- 
ing to the Boltzmann distribution is how the amount 
of fluctuations is controlled by the temperature T. At 
high temperatures, the phases will fluctuate more easily 
whereas they tend to freeze when T is lowered. In this pa- 
per, we will use units of temperature that are such that 
the critical temperature of the model in Eq. (|6]), with 
Pi^ = 1 and Ai^ = 0, is Tc — 2.2. By inspecting the en- 
ergy Eq. we see that in a Monte Carlo computation, 
the density profile Pi^ effectively works as the inverse of 
a position dependent effective temperature, such that 



T 

Teff = — , 

ill. 



(11) 



whence we expect phase fluctuations to depend strongly 
on the density profile of the condensate. This is an ex- 
act statement within the phase-only approximation. In 
particular, this means that the trapped condensates effec- 
tively are "warm" (in the sense of being close to the con- 
densation temperature) wherever the ground state den- 
sity is low. In regions of low density we therefore expect 
more violent vortex fluctuations. As will be shown below, 
this is typically the case close to the edge of the trap, but 
may also be true close to the center of the trap for anhar- 
monic trapping potentials with a dip in the condensate 
density at the center of the trap. 

We define the process of going over all sites once as one 
sweep and measure the Monte Carlo time in units of the 
sweeps. Initially at each temperature, all realizations of 
the system are thermalized with at least 100 000 sweeps 
to make sure they fluctuate around the equilibrium state 
before any measurements are made. To calculate thermal 
averages, we sample the configuration every 100th sweep. 
The numerical grid is cubic with a linear extension L = 
72 (for the case of the uniform cylinder, we also consider 
L = 36). 



A. The helicity modulus 

Phase coherence in a vortex system is probed by com- 
puting the helicity modulus T^, equivalently the su- 
perfluid density, defined as the change in free energy 
F = —T\nZ as a result of a phase twist A applied in 
the the following way. 



E[A] 



^cos{A^Oi-Ai^-A^). (12) 



1 if < 0, 

^-AE/T ifA^>0. 



The expansion of the free energy is even in A, which 
can be viewed as a change in the boundary conditions of 
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the system. For small deviations from periodic boundary 
conditions, the leading behaviour is quadratic and the 
helicity modulus is given as the coefficient to the second 
order term, 



T = 

J- //, o 



1 5^F 




• (13) 



and defining the modified helicity modulus as follows, 



Here, Yl' is over all sites where ^{viz) is nonzero (de- 
pending on Ri and R2) and TV' is the number of these 
sites. The local T z{Ri^R2) provides a measure of the 
phase coherence more locally than the global quantity, 
and will therefore be used to investigate the character of 
the vortex state at different positions in a trapped Bose- 
Einstein condensate. 



The twist can be applied in any direction, and with- 
out rotation the response in terms of is equal for all 
ji = x^y^z. When the temperature is increased from low 
to high, thermal fluctuations gradually destroy phase co- 
herence. Consequently, the renormalized superfluid den- 
sity continuously evolves from a finite value to zero at 
some critical temperature. On the other hand, in a rotat- 
ing system the helicity modulus will be different for the 
X- and ^/-directions than along the axis of rotation. Both 
cases are however important since T^^ and carries in- 
formation on numerical pinning. In the computations, we 
choose the vector potential so that V x A = (0,0,27r/), 
where / is the number of rotation-induced vortices per 
numerical grid-plaquette in the x?/-plane. For technical 
reasons, we restrict the ffiling fraction / < k/LF' with 
k = 1,2,3,..., but at the same time the density of vor- 
tex lines should not be too high in order to avoid artificial 
pinning to the underlying numerical grid^. This can be 
probed by the helicity modulus in the transverse direc- 
tions. A zero value of T^^ and means that the vortex 
line lattice is free to move translationally with respect to 
the grid. For very low temperatures, there will always be 
pinning, but the pinning should disappear well below the 
temperature at which the lattice melts. This melting on 
the other hand, is characterized by a discontinuous jump 
in the superfluid density measured parallel to the vortex 
lines, i.e. a jump in to zero. 

In non-uniform systems, we encounter a problem with 
the above definition of the helicity modulus, since Pi^ 
equal to or close to zero at the system edges will favour 
fluctuations at all temperatures in these regions. The 
global mixes together information on the amount of 
fluctuations in all regions, and the interpretation is there- 
fore less useful. To obviate this difficulty, we introduce 
a modified helicity modulus in z-direction, defined in a 
selected region between two cylinders of radii Ri and R2 . 
We do so by applying a twist 



Az if Ri < < R2, 
otherwise. 



(14) 



IV. UNIFORM SYSTEMS 

As a warmup to the results to be presented below, 
we ffist consider two cases of uniform systems, namely 
the infinite uniform system and a cylindrical confinement 
geometry with uniform density. The former in particular 
allows connections to be made to the vast literature on 
vortex physics of extreme type-II superconductors such 
as high-Tc superconductors. 



A. Infinite uniform system 

We begin with a review of the Monte Carlo results for 
a uniform system, i.e. with Pi^ = 1 in Eq. ([6]), and we 
first consider the non-rotating case. This is the standard 
uniform 3D XY model, which has been studied exten- 
sively elsewhere, see for example Ref. 34.35.36.37.38. At 
low temperatures, the phases tend to align as spins in 
a ferromagnet, and the system is frozen in a stiff state 
where any change in boundary conditions is associated 
with a large response in free energy. Consequently, the 
helicity modulus is close to unity. Equivalently, the su- 
perfluid density is close to the ground state density Pi^. 
In Fig. [1] is shown the helicity modulus along the 
z-direction, but T^^ or Ty would give similar results. As 
the temperature is increased, the relevant phase fluctu- 
ations start to appear as vortex loops, in numbers that 
gradually increase with temperature. This reduces the 
superfluid density. However, not until the vortex loops 
loose their line tension (free energy per unit length) at the 
critical temperature Tc — 2.2, will vanish completely. 
The phase transition is continuous and accompanied by 
a diverging length scale in the thermodynamic limit and 
a developing singularity in the specific heat Cy as shown 
in Fig. [TJ High precision measurements of the critical 
exponents can be found in Ref. 51 

In the second row of Fig. [U we show some snapshots 
from the computations where the vortices have been cal- 
culated via Eq. ([7j) and plotted in a 3D volume. The 
vortex radius is chosen to be 0.4 times the grid spac- 
ing for a convenient visualization, and this should not be 
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T=1:0Q- T = 1.67 T=2.50 



FIG. 1: (Color online) Monte Carlo results for a uniform sys- 
tem with periodic boundary conditions. The helicity modulus 
(superfluid density) along the 2;-direction is plotted (+), 
but Tec and Ty give equal results. The specific heat Cv/L^ 
(x) has a peak at Tc, where Tz vanishes. In the second row, 
random snapshots of the vortex configurations are shown for 
three temperatures. These are sections with 16^ lattice sites 
cut out from the 72^ system. 



associated with the true core size. Additionally, the vor- 
tices (especially the vortex lines in the rotating system we 
present below) exhibit sharp bends at short length scales. 
These bends result from the numerical grid. In fact, the 
vortices shown in the images are splines; the precise form 
of the vortices in the computations is straight line seg- 
ments connected in perpendicular corners. Nevertheless, 
the model has proved to be accurate for vortex fluctu- 
ations at scales larger than the grid spacing24i^i^i^i2^. 
Hence, the 3D snapshots provide useful hints to the state 
the system is in at different temperatures. Only occa- 
sional and small vortex loops appear when the temper- 
ature is low, but eventually they fill the whole volume. 
Above Tc, the system is in a state dominated by a tan- 
gle of tensionless vortex loops of all sizes, and there is no 
phase coherence or superfluidity. 

In a rotating system, the scenario is different. The 
ground state is the famous Abrikosov lattice^^^^^, where 
straight rotation-induced vortex lines arrange themselves 
in a triangular pattern, though with some defects due to 
the incommensurable underlying square numerical grid. 
Here, we present computation results from a system with 
filling fraction / = 1/36, and the structure of the vortex 
line lattice can be seen in Fig. [2l In the upper row is 



shown the structure function. 




where n(r^) is the local vorticity given by Eq. ([7j) and the 
sum runs over the positions of all plaquettes in the xy- 
planes. The r^- and k^-vectors are perpendicular to the 
axis of rotation. The structure function exhibits sharp 
peaks for the characteristic Bragg vectors of the vortex 
line lattice when the temperature is low and the system 
close to its ground state. As the effect of thermal fluc- 
tuations sets in and the vortex lines gets more prone to 
bending, the Bragg-peaks are weakened, and eventually 
there is a transition from a sixfold-symmetric structure 
to a ring structure in S'(k^). This is the hallmark of a 
vortex liquid phase where the vortices still posess a line 
tension which disappears in a crossover transition at an 
even higher temperature. The second row in Fig. [2] is the 
real-space equivalent of the above, where the local vortic- 
ity n{r±) is integrated along the 2:-direction so that closed 
vortex loops will be cancelled or averaged out. Both the 
Fourier- and real-space versions of the structure func- 
tion are thermal averages calculated from 100 000 Monte 
Carlo sweeps. In the left panels of the second row, bright 
spots correspond to straight lines in relatively stable po- 
sitions. Higher temperatures result in increased bending 
of the lines, and the bright spots develop into smeared- 
out regions until it is no longer possible to distinguish 
individual vortex lines in the rightmost panel. Insight 
into the bending mechanism can be obtained from the 
3D sections in the third row, which are snapshots cor- 
responding to those of the non-rotating case in Fig. [H 
In principle, such snapshots from a Monte Carlo compu- 
tation could produce any possible configuration, but a 
state far from equilibrium is highly unlikely. We there- 
fore assume the pictures to be representative and useful 
indications of the system's state at a given temperature. 
Compared to the non-rotating system, we see that there 
is much going on in terms of vortex fluctuations even at 
low temperatures like T = 1.00, in the units we defined 
below Eq. (p!Q|) . The energy cost associated with elemen- 
tary vortex excitations, i.e. a closed vortex-loop, is less 
when there are vortex lines already present. A loop and 
a line can merge to produce a bend in the vortex line 
and enough bends will result in a melting transition at 
Tm' This is succeeded by a vortex loop blowout inside 
the vortex-liquid phase similar to what happens in the 
non-rotating case. 

The qualitative picture above is supported by the qan- 
titative measurements presented in Fig. [3l The vortex 
loop blowout no longer corresponds to a phase transi- 
tion, but a remnant crossover is still indicated by a peak 
in the specific heat at a temperature T > Tm- A finite 
size scaling analysis would however reveal that there is 
no criticality associated with this broad peaki^. On the 
other hand, the melting of the vortex line lattice is a first 
order phase transition characterized by a discontinuity 
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FIG. 2: (Color online) For each temperature, given below the columns, the structure function 5'(k^) is shown in the upper 
row (5'(0) is removed, and in the two rightmost images the colorscale is magnified by a factor 25). There are sharp peaks for 
the characteristic Bragg vectors at the lowest temperatures before a transition to a ring structure corresponding to a vortex 
liquid phase at T = 1.67. In the second row we show the real-space equivalent to the structure function, the xy-positions of 
the vortices, integrated over ^-direction. The averages of the first two rows are calculated from 100 000 Monte Carlo sweeps. 
In the third row, we show sections of size 16^ from snapshots of the vortex configurations. 
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FIG. 3: (Color online) The helicity modulus Tz (+) along 
the axis of rotation is cut off and vanishes when the vortex 
lattice melts. In the transverse direction, a zero = T^, Ty 
(□) indicates that there is no pinning of the vortices to the 
numerical grid. The rounded peak in Cv/L^ (x) is a rem- 
nant of the vortex-loop blowout transition in the non-rotating 
system. 



in the helicity modulus (or superfiuid density) Tz along 
the axis of rotation in the thermodynamic limit. In a 
finite system like the one we have simulated, the drop 
to zero is continuous, but compared to a non-rotating 
system, it is much steeper and takes place at a lower 



temperature. Finally, note that the helicity modulus in 
the transverse direction is zero for temperatures well 
below the melting transition, indicating that the vortex 
line lattice is not pinned to the numerical grid. 

B. Cylindrical container 

We next consider the case of a uniform cylinder. This 
gives a ground state density profile P^^ illustrated in Fig. 
m Figs. [5] and [6] show the results from computations 
of vortex matter in a cylindrical container with such a 
density profile, given by 

P^^=e{r^^-R). (17) 

where is the Heaviside step function Q{x) = 0,x < 0; 
Q{x) = 1, X > 0. We have used two different system sizes 
L, but with the same filling fraction / = 1/36 to see 
how the radius R = L/2 — 2 affects the ordering of the 
vortices. At low temperatures, the computations repro- 
duce orderings with circular distortions of the vortex line 
lattice near the container wall, as predicted for ^He in a 
zero-temperature treatment of the problem^^. For a large 
number of vortices the system reacquires the hexagonal 
lattice symmetry away from the wall, see Fig. [5] (bottom 
row). Increasing the temperature in the case of small 
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FIG. 4: (Color online) Radial density profile Pi/j^ for the case of 
a uniform sylinder. The density is uniform in the ^-direction. 
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FIG. 5: (Color online) xy positions of vortices in a cylindrical 
container integrated over ;2-direction, and averaged over every 
tenth of a total of 5 • 10^ MC sweeps. Top and bottom rows 
have L = 36 and L = 72, respectively. At T = 0.5, 1.0, we 
discern circular ordering close to the cylinder wall combined 
with a hexagonally ordered state closer to the center. At 
T = 1.25, 1.67 we observe dominance of angular fluctuations 
closest to the edge. 



number of vortices (top row of Fig. [5|) the dominant 
vortex fluctuations are associated with angular displace- 
ments, while radially the vortex density remains ordered. 
For the largest system, with many vortices, (bottom row 
of Fig. [5]) we find dominance of angular fluctuations only 
for the vortices situated close to container wall, while the 
center of the system does not display this phenomenon. 

The crossover to a uniformly molten vortex system oc- 
curs in both cases only at a higher temperature. The 
two-step thermal crossover in the vortex pattern we find 
is analogous to that in two dimensions where vortices are 
point-like objects (see e.g. Ref. [HH). There is, however, 
a principal difference in our case, since in three dimen- 
sions the vortex line lattice melting is accompanied by 
significant vortex bending fluctuations. 

Inspecting the modified helicity modulus T z{Ri^ R2) 
(Fig. [6|) for different regions (i?i,i?2) inside the cylin- 
der, and comparing these results to those of a uniform 




T= 0.50 



FIG. 6: (Color online) Results for Tz{Ri, R2) in a cylindrical 
container. In panels (b)-(e) are shown the helicity moduli 
fz{Ri,R2) for radii Ri,R2 = 0, i?/4, 2i?/4, 3i?/4, i?, as in- 
dicated by the white circles in the images on the left, and 
compared to the results for the extended uniform system, 
panel (a). The upper curve (+) is the helicity modulus with- 
out rotation, while the lower curve (x) is the results with 
rotation-induced vortices present (filling fraction / = 1/36). 
All regions in the cylindrical container behave almost as the 
uniform system. 



system (panel (a)), we find only small differences. The 
cylindrical system behaves as the uniform one, and the 
circular distortions do not seem to affect the superfluid 
density. One could argue that the drop to zero is a little 
more rounded in the rotating non-uniform system, but 
this can be explained by the fact that T z{Ri-,R2) is cal- 
culated for smaller subsystems and that the finite size 
effects necessarily should be more severe here. In fact 
the cylindrical container is just a uniform system with 
different boundary conditions, and in a sufficiently large 
system the boundary effects are irrelevant. 
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V. NON-UNIFORM SYSTEMS 

In this section we present Monte Carlo results from 
systems with non-uniform density profiles Pi^. The fi- 
nite size of the systems is now a wanted feature and 
closer to the real situation with ultra-cold atoms, and we 
do not need to do finite size scaling. Strictly speaking, 
the only possibility for transitions are of crossover na- 
ture. To reduce the surface effects and because we model 
elongated systems, we employ periodic boundary condi- 
tions in the z-direction. Global quantities such as the 
ordinary helicity modulus and the structure function 
5'(k^) have no rigorous meaning in these systems, and we 
rather use local versions like the modified helicity mod- 
ulus T z{Ri, R2) and the local vorticity n{r±). The 3D 
snapshots are also useful indicators on the mechanisms 
involved. 



A. Harmonic trap 

For a system in a harmonic trap, we choose 
a density profile according to the Thomas-Fermi 
approximatio n^^i^^i^^ with the shape of an inverse 
parabola, Pi^ = (^) ©(^ ~ H/i), where 9 is the 
Heaviside step function Q{x) = 0, x < 0, 0(x) = 1, x > 0, 
and 

P^{x) = l-x\ (18) 

The density gradient in a trap can alternatively be viewed 
as an effective temperature gradient in a uniform system, 
see Eq. [TTJ It is clear that for low, but finite, actual tem- 
peratures T, there will be a finite area near the edge of 
the cloud which effectively would be at a high enough 
temperature to feature an annulus of tension-less tangle 
of vortices. This is a phase where the vortex line tension 
has vanished through the proliferation of vortex loops. 
The ground state density profile P^^ is shown for the har- 
monic trap in Fig. [71 The question is whether there is 




FIG. 7: (Color online) Radial density profile Pi^ for the case 
of a harmonic trap. 

also a vortex liquid region in between the ordinary vortex 



line lattice and the tensionless annulus, which represents 
the true boundary of the condensate. The 3D snapshots 
in Fig. [8] illustrate how the vortices are stiffer in the cen- 
tral part than further out towards the edge where there is 
an annulus with a tangle of tension-less vortices. At the 
higher temperature on the right, this region has grown, 
but simultaneously the amount of bends in the vortex 
lines in the center has increased. 

Further insight into the stability of the vortex line lat- 
tice can be obtained from the z-integrated vorticity in 
Fig. [9l The top row consists of snapshots and already 
here it is easy to separate the ordered lattice region from 
the disordered one, since straight vortex lines are seen as 
bright spots while bent vortices result in smeared spots or 
even smeared regions. This observation may be related 
to experiments, where at least for non-equilibrated vor- 
tex systems the ^-integration renders vortices essentially 
indistinguishable^^'^^. We have used a filling fraction 
/ = 1/36. Taking parameters from Ref. 58, and using 
ft = {h/M)Ny/2iTR^ where Ny is the number of vortices 
in the trap, we find Vt ~ lOOHz. Since uo^ ~ SOOHa^, 
this puts us well outside the LLL regime. 

We find that a well ordered vortex line lattice extends 
over most of the system at T = 0.50. Note also the ab- 
sence of circular distortions for the vortices at the edge of 
the system, as opposed to the situation in the cylindrical 
container. At T = 1.67 we can still distinguish 2 — 3 cen- 
tral vortices where the density is the highest in the snap- 
shot. However, in the thermal averages of the second and 
third row these vortices are no longer possible to detect 
due to their thermally fluctuating positions. The ther- 
mal averages are created by averaging snapshots like the 
ones in the first row over 100 000 Monte Carlo sweeps in 
the second row and 500 000 sweeps in the third. Indeed, 
in a finite system, the averaging will eventually produce 
a complete smearing even in the center of the trap since 
there is only a finite energy barrier to translate or rotate 
even a perfect vortex line lattice. Signatures of this ef- 
fect can be seen in the difference between the T = 1.25 
pictures of the second and third row of Fig. [9l However, 
it is important to note that the Monte Carlo time scale 
of the fluctuations we observe in the ordered regions, is 
dramatically larger than those related to the fluctuations 
in the disordered ones. Therefore, it does make sense to 
speak of ordered and disordered regions in these systems. 

We then again turn to the question of the character of 
the vortex state in the disordered region and the possi- 
bility of a vortex liquid layer here. For this we use the 
modified helicity modulus Tz{Ri^R2) and compare the 
results to the extended uniform system as we did for the 
cylindrical container. The results are shown in Fig. [TOt 
where both rotating and non-rotating systems are con- 
sidered. If one were to observe no appreciable difference 
in the temperature dependence of the helicity with and 
without rotation, one would conclude that the demarca- 
tion line seen in the images of Fig. [9] separates an ordered 
region from tension-less vortex tangle, with no discernible 
vortex liquid region. 
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FIG. 8: (Color online) Snapshots of vortex configurations in a rotating trapped BEC at T = 0.5 (left figures) and at T = 1.0 
(right figures). The top row shows a selection of 16 out of 72 layers in z direction. The bottom row shows smaller selections in 
the xy plane, but 32 out of 72 layers in z direction. Fluctuations are minimal in the trap center, and increase towards the edge 
of the trap. A distinct front separating regions of ordered and disordered vortices is easily identified. 
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FIG. 9: (Color online) x2/-positions of vortices in a trapped 
BEC integrated over 2;-direction. In the top row, we show 
snapshots, while the middle and bottom rows show averages 
of 10^ and 5 • 10^ MC sweeps, respectively. Every tenth con- 
figuration has been sampled. This provides information on 
the stability of the ordered region and the evolution of the 
disordered region as T varies. 



The measurements indeed show that the presence of a 
rotation significantly reduces the temperature at which 
Tz vanishes. This reduction relative to the case of no 
rotation decreases with increasing Ri^ R2. That is, panel 
(d) is similar to panel (c) (no trap), whereas in panel 
(g) there is little difference between Tz{Ri, R2) with and 
without rotation. Thus, for the latter case the presence 



of vortices essentially does not influence T^, and the de- 
struction of superfluid density is driven by the prolifera- 
tion of vortex loops. In panel (a) and (b) we have plotted 
T z as a function of the distance = {Ri + i?2)/2 from 
the center of the trap for different temperatures. The 
same conclusions may be drawn here. For the largest 
r±_ there are small differences between (a) and (b), but 
closer to the center, the helicity modulus suddenly drops 
to zero when the temperature is increased in the sys- 
tem with rotation (b). This corresponds to the melting 
transition. The reason why we only find such character- 
istics in the central part, is that the density is almost 
uniform here. Further out, the density gradient is sim- 
ply too large and the vortex line lattice crosses directly 
over to the tension-less tangle with no visible tension-full 
vortex liquid region in between. This is consistent with 
the experiments showing a very regular edge structure for 
systems with large number of vortice a^^i^^i^^i^^i^^i^'^i^^ . 



B. Anharmonic trap 

Experimentally, the traps that are used to con- 
fine the Bose-Einstein condensates may also be made 
anharmonic^^. To study the effect of a harmonic 
plus quartic trapping potential, we will use a modified 
Thomas-Fermi density distribution which varies in the 
x?/-plane as ai + 0^2 + 0^3 In an experiment 
with a rotating atomic gas, the ratio a^/a2 depends 
on the rotation frequency, but for technical reasons we 
cannot change Vt during the computations. Thus, we 
choose a fixed number of vortices by specifying / and 
vary only the temperature during computation runs for 



11 



No rotfltion 



With rotation 




0.2 0.4 0.6 0.3 




0.2 0.4 0.6 O.S 



-1 1— 



H =^ 



4- 



.0.5 



1.5 



2.5 



fixed as/a2-Tditios. To be precise, we iiave ciiosen tiie fol- 
lowing density distribution Pi^ ~ ©(^ ~ ^i/i)' 
where 6 is the Heaviside step function Q{x) = 0,x < 0, 
Q{x) = 1, X > 0, and 



^^^^""^ {l^ax^-{l^a)x^}. (19) 



4(1 + a) 



This ground state density profile has the property that 
its maximum is always unity. Consequently, the regions 
with the highest density should be comparable to both 
the uniform system and to the center of the harmonically 
trapped system at any given temperature. The value of 
a determines the shape of the density profile. Note that 
a = corresponds to a pure quartic trapping potential. 
We present results for this in the upper row of Fig. [12] 
for filling fractions / = 1/36 (left) and / = 1/18 (right). 
Such systems have a large density gradient close to the 
edge, but is on the other hand more uniform in the in- 
ner parts than the harmonic case. The effect is a sce- 
nario which fits in between the cylindrical case and the 
harmonically trapped system. There are slight angular 
distortions of the outermost vortices for / = 1/36, and 
an increased possibility for vortex lattice melting in the 
center as in an extended uniform system. 

The generic ground state density profile P^^ for the an- 
harmonic trap is shown in Fig. [TTl where the parameter 
a in Eq. ([19]) is taken to be a = 2, such that the density 
profile has a distinct dip in the trap centre. 



T= 0.50 




FIG. 10: (Color online) Results for fz{Ri,R2). The two top 
panels show computation results for the thermal depletion of 
the superfluid density in a harmonically trapped condensate 
{r± is the distance from the center of the trap) at the temper- 
atures T = 2.50 (the lowermost curve), T = 2.00, T = 1.67, 
T = 1.25, T = 1.00, T = 0.50. The uppermost curve is 
the pure ground state density profile Pi^. In panels (c)-(g), 
the upper curve (+) is the helicity modulus without rotation, 
while the lower curve ( x ) the helicity modulus with rotation- 
induced vortices with filling fraction / = 1/36 as functions 
of temperature. Panel (c) shows Tz for a cubic uniform sys- 
tem with periodic boundary conditions. The remaining pan- 
els (d)-(g) show t40,i?/4), fz{R/4,2R/4), t,(2i?/4, 3i?/4), 
and Tz{3R/4, R), respectively. The vortex plots on the left 
(obtained at T = 0.50) defines the radii Ri and R2 as white 
circles. 



FIG. 11: (Color online) Radial density profile Pi^ for the case 
of an anharmonic trap with a local minimum of the ground 
state condensate density at the center of the trap. 



In the second row of Fig. [121 we have chosen a = 2 
so that the density profile has a local minimum in the 
center. For low temperatures, the vortex configuration is 
close to that of the harmonic system, see Fig. [9] Then we 
notice a decreased visibility of the central vortices from 
around T = 1.25, a feature encountered in an experi- 
ment reported by Bretin and coworkers in 2004^^, where 
they increased the rotation frequency above the trap fre- 
quency. They speculated that this could be explained by 



12 



•.v/.\v. 

• • • • « 
, k » • • 

• • • • 


2^^.^.1.11 

/ • • • • * * 

f r* • * < ■ • 


T-1.43 T = 1.67 

:v:v.v. ;v::;}: 


T,~i--.00 


T =,1.11 T = 1.25 T = 1.33 T = 1.43 


T.W:iJ)0 

/ • • * • . 

(••:'••.••::) 




T .^-1.25 
••'•*•". 
'-'1'. *,"•" • 


r = i.67 








T = 1.33 
^^^^ 


T = 1.43 



/ = 1/36 / = 1/18 



FIG. 12: (Color online) x^-positions of the vortices in a pure quartic trapping potential (top row) with filling fractions / = 1/36 
(left) and / = 1/18 (right). In the second row, the vortices are trapped by a harmonic plus quartic potential, with a = 2 in 
Eq. ()19|) . Due to lower density, the visibility of the vortices close to the rotation axis is reduced for the highest temperatures. 
All images are averages over 100 000 Monte Carlo sweeps. 



bending effects of the vortices. Our computations sup- 
port this view. 




T=1.25 T 

FIG. 13: (Color online) Results for Tz{Ri,R2), comparing 
the low- and high-density regions in the top and bottom pan- 
els respectively. The lower curve (x) in both panels, as well as 
the images on the left, are calculated for a system with filling 
fraction / = 1/18 and anharmonic density profile defined by 
a = 2 (see Eq. (p^ ). Also included, are results for the same 
density profile, but with no rotation-induced vortices (upper 
curve, (+)). 



Fig. [13] gives examples of the modified helicity modu- 
lus Tz{Ri^ R2) for two distinct regions of a system with 
/ = 1/18 and a = 2\ One covering the region closest to 
the rotation axis, and one around the peak in the den- 
sity. We have also included results for a non-rotating 
system, but with the same density profile. The vortex 
lines clearly forces T z{Ri^R2) to vanish in both regions 
for lower temperatures than without rotation, indicating 
a melting of the vortex line lattice. Additionally, we see 
that the region with the highest density remains in the 
vortex solid state up to higher T than the inner part. 



The central vortex line lattice melts easier because in 
this region there is a larger tendency for the vortex lines 
to have bending fluctuations due to the higher effective 
temperature, in the sense of Eq. (pTj) . 



VI. CONCLUSIONS 

In this paper, we have considered thermal fluctuations 
of vortex matter in cylindrical traps with varying non- 
uniform model ground state densities of the condensate. 
As a benchmark we have also performed Monte Carlo 
computations of a uniform system and a cylindrical ge- 
ometry with uniform ground state condensate density. 
We have in all cases taken care in ensuring that the he- 
licity moduli in the direction transverse to the direction 
of rotation have vanished at temperatures well below the 
temperatures where the helicity moduli along the direc- 
tion of rotation vanish. The purpose of this is to mimick, 
as well as possible, phase transitions in the BEC-vortex 
system in a continuum and to eliminate the artificial pin- 
ning effects that are inevitably present due to the numer- 
ical grid. Note however, that the model defined on a nu- 
merical grid could have a physical realization in terms of 
rotating Bose-Einstein condensates on an optical lattice. 
In such a system there is a real frustration of the vor- 
tex system due to commensuration effects. While this 
presents a very interesting problem in its own right, a 
study of such effects has not been the purpose of this 
paper. 

Uniform system: In a uniform system, the 3D XY spe- 
cific heat anomaly which the system exhibits in the ab- 
sence of rotation, is rounded to a broad non-critical peak 
when rotation is present. This is the superfluid/BEC 
counterpart of the well-known finite-field remnant of the 
zero-magnetic field superconductor-normal metal transi- 
tion when a supercondutor is subjected to a finite mag- 
netic field. In principle, it is also possible to extract a 
specific heat anomaly associated with the melting of the 
BEC-vortex lattice, although this has not been done in 
the present paper. The vortex fluctuations of importance 
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in this case are transverse phase fluctuations of the su- 
perfluid order parameter. 

Uniform cylindrical confinement: In a cyhndrical con- 
flnement with uniform ground state density proflle, we 
have investigated the thermal fluctuations of the ordered 
vortex lattice in the conflnement geometry as the tem- 
perature is increased. We flnd that the vortex lattice 
is flrst excited close to the circular edge of the conflne- 
ment, and that the fluctuations in the vortex positions 
are directed along the perimeter of the cylinder. As the 
temperature is increased further, these ring-like thermal 
excitations creep inwards towards the center of the cylin- 
der until the vortex lattice also at the center crosses over 
to a liquid phase. It is important that the inhomogeneity 
observed in the vortex fluctuations in this case is not due 
to the fact that the condensate effectively is warmer (in 
the sense of Eq. (pT]) ) at the edge than at the center. 
The effective temperature is uniform throughout the sys- 
tem. The reason for the observed inhomogeneity of the 
fluctuations is that the local environment around each 
vortex is different at the edge than at the center. In par- 
ticular, the vortex- vortex correlation effects that impede 
thermally assisted vortex motion are smaller at the edge 
than at the center. The coordination number of each vor- 
tex in the vortex lattice is larger at the center than at the 
edge. The corresponding interaction energies, and hence 
Coulomb barrier, that must be overcome to produce vor- 
tex motion, is therefore smaller at the edge of the con- 
flnement than at the center. Moreover, a transverse (an- 
gular) motion-pattern of the vortices are preferred com- 
pared to a radial motion-pattern, for the following rea- 
son. In a superfluid/BEC, the vortex interactions are un- 
screened (anti) Biot-Savart interactions. Therefore, the 
vortex ensemble constitutes an incompressible system in 
a superfluid/BEC. Hence, vortex motion is a collective 
process involving many vortices that have to be 'pushed' 
out of the way in order to pave the way for one vortex. 
Motion along the perimeter may be realized by moving 
all vortices in the outermost ring collectively, while ra- 
dial motion involves a rearrangement of all vortices in 
the system. The former is obviously a collective process 
that is easier to accomplish than the latter. As a corol- 
lary, we may infer that this ring-like excitation pattern 
in a cylindrical conflnement may change in the case of a 
moderate type-II superconductor, where effective screen- 
ing of vortex interactions is a consequence of a fluctuating 
gauge fleld. This results in a compressible vortex system 
which to a much larger extent will allow radial as well as 
angular fluctuations in vortex positions. 

Harmonic and anharmonic traps: For a Bose-Einstein 
condensate in a magnetic trap, we have focused on two 
types of trapping potentials, namely i) a harmonic trap 
giving a Thomas- Fermi ground state density proflle of the 
BEC, and ii) an anharmonic trap in which we could vary 
the relative weights between a quartic and a quadratic 
term. The latter trap naturally leads to a modiflcation of 
the inverse-parabolic Thomas- Fermi density proflle, and 
in particular it is possible to induce a ground state super- 



fluid density with a local minimum at the center of the 
trap. The thermally driven vortex-excitations in these 
systems are fundamentally different from the case of a 
uniform cylinder. The reason is that the vortex matter in 
the trapped BECs effectively have a highly non-uniform 
temperature, in the sense of Eq. (pT]) . This is due to the 
rapidly decreasing density proflle close to the edge, with 
a maximum effective temperature gradient at the edge of 
the trap. Hence, this promotes vortex-loop excitations 
as the edge of the trap is approached, thus inducing a 
vortex-matter phase where the line tension (free energy 
per unit length) of the vortex lines has vanished and the 
vortex lines effectively have lost their directionality in 
the direction of the rotation vector of the system. This 
gives a distinct region at the outer edge of the trap where 
vortex-loop induced violent vortex-line fluctuations wash 
out the image of each individual rotation-induced vortex 
line. This contrasts sharply with the images of ring-like 
collective directed vortex-line excitations in the uniform 
cylinder. The main difference between the harmonic and 
anharmonic trap is that for a harmonic trap the effec- 
tive temperature of the system decreases monotonously 
towards the center of the trap, such that a monotonous 
evolution of the system from a vortex tangle to a vor- 
tex lattice at the center is observed when heated. For 
an anharmonic trap with a dip in the ground state con- 
densate density at the center of the trap, we effectively 
have a local increase in the temperature of the system as 
the center of the trap is approached. Hence, we can have 
a tensionless vortex tangle at the center of the trap as 
well as at its edge. Our computations therefore provide 
a flnite-temperature extension of the zero-temperature 
phase diagram for anharmonic traps presented in Ref. 
[59I . According to Refs. [6Q||61| , the experiment by Bretin 
et al.-^^ takes place in the same regime, in which the den- 
sity and interaction strength are too high to obtain a 
giant vortex with multiple quantization. We believe the 
experiment can be described by our model. 



To summarize, we studied effects of density inhomo- 
geneity and flnite-size in a model system describing a 
trapped Bose-Einstein condensate. Although flniteness 
of the system and density inhomogeneity prohibit use 
of the notion of a true phase transition between differ- 
ent aggregate states of vortex matter (known to occur in 
three dimensional extended model systems), we flnd that 
nonetheless one can have various quasi-states of vortex 
matter surviving in a rather robust form at flnite length 
scales in traps. Recent progress in Bose-Einstein con- 
densates has resulted in the availability of systems with 
various symmetries and multiple components where one 
can anticipate new states of vortex matter. This study 
suggests that predictions of a theory based on a uniform 
density might nonetheless have rather robust flnite-size 
realizations in actual inhomogeneous trapped systems. 



14 



Acknowledgments 

This work was supported by the Research Coun- 
cil of Norway, Grant Nos. 158518/431, 158547/431 
(NANOMAT), 167498/V30 (STORFORSK), the Na- 



tional Science Foundation, Grant No. DMR-0302347, 
Nordforsk Network on Low-Dimensional Physics. SK and 
EB acknowledges the hospitality of the Center for Ad- 
vanced Study at the Norwegian Academy of Science and 
Letters, where part of this work was done. 



^ I. R. Coddington, Vortices in a Highly Condensed Bose 
Gas, PhD Thesis, University of Colorado, (2004). 

^ M. Greiner, Ultracold quantum gases in three-dimensional 
lattice potentials, PhD thesis, Ludwig Maximilians Univer- 
sitat Muenchen, Germany (2003). 

^ C. Regal, Experimental realizations of BCS-BEC crossover 
physics with a Fermi gas of atoms, PhD thesis, University 
of Colorado (2005). 

^ M. W. Zwierlein, High- Temperature Superfluidity in an Ul- 
tracold Fermi Gas, PhD Thesis, Massachusetts Institute of 
Technology (2006). 

^ J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 
(1973); J. M. Kosterlitz, J. Phys. C 7, 1084 (1974); ibid, 
10, 3753 (1977). 

^ For a thorough and comprehensive statistical mechani- 
cal treatment of the metal- insulator transition in the 2D 
Coulomb gas, with particular emphasis on the difficult 
problem of treating the screening properties in the criti- 
cal region, see J. S. H0ye and K. Olaussen, Physica 104A, 
447 (1980); ibid, 107A, 241 (1981). 

^ H. Kleinert, Lett. Nuovo Cimento, 35, 409 (1982). 

^ Z. Tesanovic, Phys. Rev. B 51, 16204 (1995); ibid, B 59, 
6449 (1999). 

^ A. K. Nguyen and A. Sudb0, Phys. Rev. B 57, 3123 (1998); 

ibid, B 58, 2802 (1998). 
^° A. K. Nguyen and A. Sudb0, Phys. Rev. B 60, 15307 
(1999); Europhys. Lett. 46, 780 (1999). See also A. K. 
Nguyen, R. E. Hetzel, and A. Sudb0, Phys. Rev. Lett., 77, 
1592 (1996). 

H. Kleinert, Gauge Fields in Gondensed Matter Physics, 
Vol. 1, World Scientific, Singapore 1989. 
"^^ K. Fossheim and A. Sudb0, Superconductivity: Physics and 
Applications, John Wiley & Sons, Ltd (2004), see Chapters 
9,10. 

The 1949 discussion remark of L. Onsager is well worth 
quoting here: "Finally, we can have vortex rings in the 
liquid and the thermal excitation of He II, apart from 
phonons, is presumably due to vortex rings of molecular 
size. As a possible interpretation of the A-transition, we 
can understand that when the concentration of vortices 
reaches the point where they form a connected vortex tan- 
gle throughout the liquid, then the liquid becomes nor- 
mal." L. Onsager, Nuovo Cimento Suppl., 249 , (1949). 
This idea was first put on a quantitative footing in 2D in 
Refs. i, and in 3D in Refs. 7,8,9,10. In particular, see Fig. 
3 of A. K. Nguyen and A. Sudb0, Phys. Rev B 60, 15304 

(1999) . 

^'^ J. Hove, S. Mo, and A. Sudb0, Phys. Rev. Lett., 85, 2368 

(2000) . 

P. W. Anderson, Basic Notions in Gondensed Matter 
Physics, Benjamin Cummings, London (1984). 
W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin, and 
D. E. Pritchard, Phys. Rev. Lett., 73, 2253 (1995); K. B. 
Davis, M. -O. Mewes, M. R. Andrews, N. J. van Druten, 



D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. 
Lett. 75, 3969 (1995) 

W. Petrich, M. H. Anderson, J. R. Ensher, and E. A. Cor- 
nell, Phys. Rev. Lett., 74, 3352 (1995); M. H. Anderson, 
J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. 
Cornell, Science, 269, 198 (1995). 

C. C. Bradley, C. A. Sackett, and R. G. Hulet, Phys. Rev. 
Lett., 78, 985 (1997). 

M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, 
C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 83, 
2498 (1999). 

2° C. Raman, M. Kohl, R. Onofrio, D. S. Durfee, C. E. Kuk- 
lewicz, Z. Hadzibabic, and W. Ketterle, Phys. Rev. Lett. 
83, 2502 (1999). 

2^ M. Snoek and H. T. C. Stoof, Phys. Rev. Lett. 96, 230402 
(2006). 

^2 M. Snoek and H. T. C. Stoof, Phys. Rev. A 74, 033615 
(2006). 

M. M. Parish, S. K. Baur, E. J. Mueller, and D. A. Huse, 
arXiv: 709. 1120 . (2007). 

E. Babaev, A. Sudb0, and N. W. Ashcroft, Nature, 431, 
666 (2004). 

S. Kragset, E. Babaev, and A. Sudb0, Phys. Rev. Lett., 
97, 170403 (2006). 

A. J. Leggett, Rev. Mod. Phys., 73, 307 (2001), and refer- 
ences therein. 
'^^ E. H. Brandt, Prog. Phys. 58, 1465 (1995). 

M. Greiner, C. A. Regal, and D. S. Jin, Nature 426, 537 
(2003). 

2^ C. A. Regal, M. Greiner, D. S. Jin, Phys. Rev. Lett. 92, 
040403, (2004). 

2° M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. 

Ketterle, Nature, 442, 54 (2006). 
2^ Y. Shin, M. W. Zwierlein, C. H. Schunck, A. Schirotzek, 

and W. Ketterle, Phys. Rev. Lett. 97, 030401 (2006). 
22 C. H. Schunck, M. W. Zwierlein, A. Schirotzek, and W. 

Ketterle, Phys. Rev. Lett. 98, 050404 (2007). 
22 C.H. Schunck, Y. Shin, A. Schirotzek, M.W. Zwierlein, and 

W. Ketterle, Science, 316, 867 (2007). 
2^ R. E. Hetzel, A. Sudb0, and D. A. Huse, Phys. Rev. Lett. 

69, 518 (1992). 

25 T. Chen and S. Teitel, Phys. Rev. B 55, 15197 (1997). 
2^ X. Hu, S. Miyashita, and M. Tachiki, Phys. Rev. Lett. 79, 
3498 (1997). 

2'^ P.Olsson and S. Teitel, Phys. Rev. Lett. 80, 1964 (1998). 
2^ S. Ryu and D. Stroud, Phys. Rev. B 57, 14476 (1998). 
29 G. Baym and C. J. Pethick, Phys. Rev. Lett., 76, 6 (1996). 

G. Watanabe, G. Baym, and C. J. Pethick, Phys. Rev. 

Lett., 93, 190401 (2004). 

N. R. Cooper, S. Komineas, and N. Read, Phys. Rev. A 

70, 033604 (2004). 

^2 K. W. Madison, F. Chevy W. Wohlleben, and J. Dalibard, 

Phys. Rev. Lett., 84, 806 (2000). 
^2 L Coddington, P. Engels, V. Schweikhard, and E. A. Cor- 



15 



nell, Phys. Rev. Lett., 91, 100402 (2003). 

V. Schweikhard, 1. Coddington, P. Engels, V. P. Mogen- 

dorff, and E. A. Cornell, Phys. Rev. Lett., 92, 040404 

(2004). 

N. L. Smith, W. H. Heathcote, J. M. Kmeger, and C. J. 
Foot, Phys. Rev. Lett., 93, 080406 (2004). 
I. Coddington, P. C. Haljan, P. Engels, V. Schweikhard, S. 
Tung, and E. A. Cornell, Phys. Rev. A 70, 063607 (2004). 
S. R. Muniz, D. S. Naik, and C. Raman, Phys. Rev. A 73, 
041605 (2006). 

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. 
H. Teller, and E. Teller, J. Chem. Phys., 21, 1087 (1953). 
W. K. Hastings, Biometrika, 57, 97 (1970). 
The resulting phase where one has zero helicity moduli in 
the direction transverse to the rotation vector, concomi- 
tant with a finite helicity modulus parallell to the rotaton 
vector, may be dubbed a "floating solid phase", since the 
vortex lattice "floats" as a solid while being thermally de- 
pinned from the numerical grid. Going to low enough filling 
fraction to avoid artificial pinning of the vortex lattice to 
the numerical grid, is a minimum requirement for studying 
phase transitions in vortex lattices in continuum systems. 
There are, however, interesting physical situations where 
one really has a rotating BEC defined on a discrete lattice 
system, most notably an optical lattices. (See for instance, 
H. Pu, L. O. Baksmaty, and N. P. Bigelow, Phys. Rev. 
Lett., 94, 190401 (2005); S. Tung, V. Schweikhard, and E. 
A. CorneU, Phys. Rev. Lett., 97, 240402 (2006).) In this 



case, one has competing length scales in the system, with 
associated competing energy scales. While this poses an 
extremely interesting problem, both in 2D as well as in 3D 
systems, it is beyond the scope of the present work. 
M. Campostrini, M. Hasenbusch, A. Palisetto, P. Rossi, 
and E. Vicari, Phys. Rev B 63, 214503 (2001). 
For a comprehensive review on the Abrikosov vortex lat- 
tice in type-II superconductors in a magnetic field, see G. 
Blatter, M. V. Feigel'man, V. B. Geshkenbein, A. L Larkin, 
and V. M. Vinokur, Rev. Mod. Phys., 66, 1125 (1994). 
S. A. Gilford and G. Baym, Phys. Rev. A 70, 033602 
(2004). 

L. J. Campbefl and R. M. Ziff, Phys. Rev. B 20, 1886 
(1979). 

Yu. Lozovik and E. Rakoch, Phys. Rev. B 57, 1214 (1998). 
V. Bretin, S. Stock, Y. Seurin, and J. Dalibard, Phys. Rev. 
Lett. 92, 050403 (2004). 

J. R. Abo-Shaeer, C. Raman, and W. Ketterle, Phys. Rev. 
Lett. 88, 070409 (2002). 

J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket- 
terle, Science 292, 476 (2001). 

G. M. Kavoulakis and G. Baym, New Journal of Physics 
5, 51.1 (2003). 

A. D. Jackson, G. M. Kavoulakis, and E. Lundh, Phys. 
Rev. A 69, 0536619 (2004). 

A. D. Jackson and G. M. Kavoulakis, Phys. Rev. A 70, 
023601 (2004). 



